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LONG-TERM  GOALS 

The  major  objective  of  this  proposal  is  to  calculate  the  3-D,  time  dependent  radiation  field  both  within 
the  ocean  and  in  the  atmosphere  in  the  presence  of  a  stochastically  varying  interface  which  may  also  be 
perturbed  by  sea  foam,  air  bubbles,  surfactants,  rain,  etc.  This  study  will  serve  as  the  genesis  to  the 
future  evolution  of  an  inversion  algorithm  whereby  one  could  reconstruct  images  that  have  been 
distorted  by  the  interface  between  the  atmosphere  and  the  ocean  or  the  ocean  itself.  This  study  will 
rely  heavily  on  both  the  spectral  and  polarimetric  properties  of  the  radiation  field  to  deduce  both  the 
sea  state  and  the  perturbations  produced  on  it.  A  second  phase  of  this  study  will  be  to  explore  the 
asymptotic  polarized  light  field  and  to  determine  how  much  information  can  be  obtained  about  the 
IOP’s  of  the  medium  by  measuring  it.  The  third  phase  of  this  proposal  will  deal  with  the  problem  of 
improving  image  resolution  in  the  ocean  using  some  novel  polarimetric  techniques  that  we  are  just 
beginning  to  explore.  Once  these  studies  have  been  completed  using  a  passive  source,  it  will  be  rather 
straightforward  to  extend  them  to  active  sources  where  we  can  explore  the  use  of  both  photo-acoustic 
and  ultrasound-modulated  optical  tomography  to  improve  image  resolution. 

OBJECTIVES 

The  new  Navy  initiative  is  focusing  on  one  of  the  most  fonnidable  problems  in  radiative  transfer 
theory;  namely,  calculating  the  full  3D  time  dependent  radiation  field  (with  full  Mueller  matrix 
treatment)  in  a  coupled  atmosphere-ocean  system  where  the  boundary  separating  the  two  has  both 
spatial  and  temporal  dependence.  Although  a  great  deal  of  work  has  been  done  on  obtaining  power 
spectra  for  ocean  waves,  I  know  of  no  work  that  has  yielded  similar  results  for  the  radiation  field 
within  the  ocean.  It  is  clear  that  as  long  as  the  surface  has  a  significant  effect  on  the  internal  light  field, 
it  will  leave  its  signature  on  the  radiation  field  within  the  ocean  and  the  relative  strength  of  this  field 
compared  to  the  ambient  field  will  determine  the  success  or  failure  of  inversion  algorithms.  However, 
as  we  go  deeper  within  the  ocean  we  start  to  enter  a  region  called  the  asymptotic  region  where  all 
photons  lose  memory  of  their  origin  and  the  light  field  then  remains  stationary  and  becomes 
independent  of  the  azimuthal  angle.  The  depth  dependence  becomes  simply  exponential,  i.e.  L(z+h,0) 
=  L(z,  0)  exp(-kh)  where  k  is  called  the  diffusion  exponent.  It  should  be  noted  at  this  juncture  that  this 
asymptotic  light  field  is  still  polarized  which  is  why  we  used  the  bold-faced  vector  notation.  We  were 
the  first  to  compute  the  degree  of  polarization  for  this  asymptotic  light  field  for  Rayleigh  scattering  and 
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were  able  to  obtain  an  analytic  expression  for  both  the  polarized  radiation  field  and  the  diffusion 
exponent  (see  ref  1).  In  addition,  we  were  also  able  to  set  up  a  numerical  scheme  to  compute  the 
polarized  radiation  field  as  well  as  the  diffusion  exponent  for  any  single  scattering  Mueller  matrix. 

The  interesting  feature  about  the  asymptotic  light  field  is  that  it  depends  profoundly  on  both  the  single 
scattering  albedo  as  well  as  the  phase  function  of  the  medium.  We  also  found  that  substantial  errors 
will  occur  in  both  the  ordinary  radiance  and  the  diffusion  exponent  if  they  are  calculated  from  scalar 
rather  than  vector  theory 

APPROACH 

There  are  several  stages  to  our  approach  that  we  will  enumerate.  The  sine  qua  non  for  this  entire 
project  will  be  the  development  of  a  fully  time  dependent  3-D  code  capable  of  calculating  the  complete 
radiation  field,  i.e.  the  complete  Mueller  or  Green  matrix  at  any  point  within  the  atmosphere-ocean 
system.  This  of  course  implies  that  both  horizontal  as  well  as  vertical  IOP’s  must  be  accounted  for.  It 
should  also  be  noted  that  the  code  must  be  capable  of  handling  internal  sources  as  well  in  order  to 
explore  both  fluorescence  and  bioluminescence.  At  present  there  are  several  3D  codes  that  are  able  to 
compute  various  radiometric  quantities  in  inhomogeneous  media;  however,  as  far  as  we  know,  none 
exists  which  will  couple  both  atmosphere  and  ocean  with  a  time  dependent  stochastic  interface.  One 
of  the  earliest  3D  radiative  transfer  (RT)  codes  was  developed  by  Stenhohn,  et.  al2  to  model  thennal 
emission  from  spherical  and  non-spherical  dust  clouds.  It  was  based  on  an  implicit  discretization  of 
the  transfer  equation  in  Cartesian  frames.  To  our  knowledge,  the  first  3D-scalar  RT  code  using 
discrete  ordinates  was  written  by  Sanchez3  et  al.;  however,  it  did  not  make  use  of  spherical  harmonics 
and  lacked  efficiency  and  accuracy  particularly  for  small  viewing  angles.  The  addition  of  polarization 
to  the  3D  discrete  ordinates  method  was  done  by  Hafennan4  et  al.  Almost  concurrently,  a  3D-scalar 
RT  code  was  written  by  K.  F.  Evans5  which  used  both  spherical  harmonics  and  discrete  ordinates.  This 
method  uses  a  spherical  hannonic  angular  representation  to  reduce  memory  and  CPU  time  in 
computing  the  source  function  and  then  the  RT  equation  is  integrated  along  discrete  ordinates  through 
a  spatial  grid  to  model  the  radiation  streams.  We  have  already  obtained  this  code  and  will  use  it  for 
validation  of  our  3D  scalar  Monte  Carlo  code  for  both  the  atmosphere  and  ocean  components.  Several 
Monte  Carlo  codes  both  scalar  and  vector  have  been  published  for  solving  specialized  problems  in 
atmospheric  optics  usually  dealing  with  finite  clouds  ’  ’  ’  .  Without  exception,  these  codes  are  using 
quite  primitive,  also  called  “brute  force”  ,  methods.  None  of  them  will  do  what  we  are  proposing  in 
our  approach  to  the  fully  time  dependent  3D  solution  applicable  to  both  atmosphere  and  ocean.  It 
should  be  mentioned  that  we  have  already  successfully  added  to  our  Green  matrix  Monte  Carlo  code 
the  capability  to  handle  internal  sources  such  as  fluorescence,  bioluminescence  and  even  thermal 
emission. 

Once  we  have  developed  our  3D  code  to  handle  both  the  ocean  and  atmosphere  without  the  interface 
included,  we  will  then  develop  a  3D  form  of  ID  matrix  operator  theory  that  we  worked  out  in  two 
seminal  papers  published  in  Applied  Optics10'11.  The  basic  idea  of  the  method  is  that  if  one  knows  the 
reflection  and  transmission  operators  of  say  two  layers,  then  it  is  rather  straightforward  to  get  the 
reflection  and  transmission  operators  for  the  combined  layer.  With  this  method  we  can  start  from  an 
infinitesimal  layer  and  build  large  and  even  semi-infinite  layers  in  a  rapid  way,  i.e.  if  we  start  with  a 
layer  of  thickness  A  then  in  N  steps  we  can  reach  a  thickness  of  2NA.  Another  very  relevant  feature  of 
this  method  is  that  it  will  allow  us  to  add  the  interface  to  the  “bare”  ocean,  i.e.,  one  without  an 
interface,  to  get  the  combined  ocean-interface  operator  and  then  add  this  layer  to  the  atmosphere  for 
the  final  reflection  and  transmission  operators  for  the  combined  system.  The  question  immediately 
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arises  is  why  not  solve  the  entire  system  at  the  same  time?  The  answer  is  that  by  doing  it  this  way  we 
only  have  to  use  the  adding  feature  to  combine  the  time  dependent  interface  thus  avoiding  performing 
the  entire  calculation  at  each  instant  in  time.  This  method  has  also  become  known  as  the  adding¬ 
doubling  method.  These  operators  are  effectively  the  impulse  response  or  Green  matrix  for  the  upper 
and  lower  boundaries  of  the  medium.  Therefore,  if  we  know  the  external  radiance  input  into  both 
upper  and  lower  layers,  we  can  then  obtain  the  output  at  both  the  upper  and  lower  boundaries  of  the 
combined  system.  A  pictorial  description  of  the  method  is  shown  in  Fig.  1 .  It  should  be  emphasized 
that  this  method  will  also  handle  internal  sources  as  well  such  as  bioluminescence,  fluorescence,  and 
even  thermal  emission.  This  method  can  also  handle  detectors  at  any  interior  point  in  the  medium. 
Another  bonus  of  matrix  operator  theory  is  that  one  can  easily  obtain  the  path  radiance  between  source 
and  detector  which  is  a  sine  qua  non  for  image  analysis. 

In  order  to  add  a  interface  which  is  spatially  inhomogeneous  in  the  y  direction  but  homogeneous  in  the 
x  direction,  we  will  need  the  reflection  (R)  and  transmission  (T)  operators  for  both  the  atmosphere  and 
ocean  now  as  a  function  of  time  t  and  both  z  and  y;  namely  R(y-yo,z,t,9,(j))  and  similarly  for  T.  It  is 
important  to  note  that  we  only  have  to  obtain  the  response  of  the  atmosphere  or  ocean  to  a  single  line 
source  at  the  point  yo  and  then  using  the  translational  invariance  of  the  medium  in  the  y-direction  will 
have  the  reflection  and  transmission  operators  at  every  point  in  the  y-direction.  The  only  method  we 
know  to  create  these  2-D  operators  is  Monte  Carlo.  Once  these  are  obtained,  we  can  use  the  output 
from  each  layer  as  the  input  to  the  surface  boundary  whose  reflection  and  transmission  properties  are 
either  known  or  calculable.  For  instance  if  the  surface  consists  of  just  capillary  or  gravity  waves,  then 
we  just  need  the  Fresnel  coefficients  to  give  us  the  requisite  reflection  and  transmission  operators  for 
the  interface.  Now  once  these  operators  are  obtained  then  we  can  use  extended  matrix  operator  theory 
to  get  the  final  time  dependent  radiation  field  that  a  detector  will  see.  Let  us  consider  the  simplest  case 
where  the  surface  is  1-D  and  we  know  its  power  spectrum.  It  should  be  emphasized  at  this  point  that  it 
is  not  sufficient  to  have  just  a  wave-slope  distribution  since  it  will  only  give  us  statistically  averaged 
results  for  the  radiance  field.  The  introduction  of  the  spatial  and  temporal  dependent  interface  destroys 
the  symmetry  and  makes  all  1-D  codes  essentially  useless  in  this  domain.  At  each  instant  in  time,  the 
surface  will  have  a  distinct  shape  that  will  evolve  in  time.  We  have  developed  a  method  using  linear 
filter  theory  whereby  we  can  take  an  ocean  power  spectrum  and  using  a  random  number  generator 
crate  a  realistic  surface  that  will  match  the  original  power  spectrum  and  will  still  exhibit  both 
stationarity  and  ergodicity.  Now  the  nice  feature  about  what  we  are  proposing  is  that  we  can  now 
concentrate  on  just  the  effects  of  the  surface  on  the  detectors  since  as  the  surface  evolves  in  time  so  too 
will  the  radiance  field  as  recorded  by  the  detectors.  Now  both  the  spatial  and  temporal  profiles  will  be 
constantly  changing;  however,  we  will  have  created  them  from  a  medium  which  has  been  assumed 
stationary  and  only  the  interface  produces  the  time  dependence  and  the  horizontal  spatial 
inhomogeneities,  i.e.  the  R  and  T  operators  for  both  the  atmosphere  and  ocean  need  only  to  be 
computed  once.  This  is  clearly  a  first-order  solution  to  the  more  complex  problem;  however,  it  should 
tell  us  a  great  deal  about  future  complexities  of  inversion  and  also  the  efficacy  of  pursuing  the  next 
level  of  difficulty.  If  the  surface  is  perturbed  by  foam,  bubbles,  etc.  then  these  can  be  added  and  the 
matrix  operator  theory  will  be  used  to  calculate  the  effective  reflection  and  transmission  operators  of 
the  perturbed  surface.  It  should  also  be  stated  that  this  project  is  enonnously  computationally 
intensive;  however,  the  type  of  codes  we  will  produce  are  ideally  suited  for  large-scale  parallel 
processors,  which  we  do  have  access  to. 

The  next  level  of  difficulty  is  where  we  will  use  Monte  Carlo  methods  to  compute  a  full  3-D 
distribution  of  the  time  dependent  radiation  field,  which  now  may  include  3-D  inhomogeneities  in  both 
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the  ocean  and  atmosphere.  This  will  be  a  computational  tour  de  force  requiring  a  major  new  computer 
program  that  must  be  capable  of  placing  IOP’s  of  both  the  atmosphere  and  ocean  at  each  point  in  a 
large  3-D  grid.  Matrix  operator  theory  will  again  be  used  but  it  will  now  be  much  more  complicated 
since  our  reflection  and  transmission  operators  now  become  functions  of  three  spatial  variables.  In 
fact,  the  complete  solution  to  this  problem  could  approach  the  complexity  of  the  general  circulation 
models  used  in  weather  forecasting.  Due  to  the  large  volume  of  data  that  will  be  generated,  we  will 
clearly  have  to  develop  methods  to  easily  display  animated  sequences  of  this  time  dependence. 

There  projects  were  worked  on  by  Dr.  Pengwang  Zhai,  who  left  for  NASA  Langley  Research  Center  in 
August  of  2008,  and  are  now  being  worked  on  by  Dr.  Yu  You,  a  postdoctoral  research  associate. 
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Interface  maps  downward  radiance  from  atmosphere  into  both  an  upward  and 
downward  radiance  then  used  as  input  into  both  atmosphere  and  ocean 


Interface  maps  upward  radiance  from  ocean  into  both  an  upward  and 
downward  radiance  then  used  as  input  into  both  ocean  and  atmosphere 
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Fig.  1.  Schematic  representation  on  the  use  of  matrix  operator  theory  to  calculate  the  time 

dependent  radiation  field  within  the  ocean 


WORK  COMPLETED 

a)  We  have  applied  the  hybrid  matrix  operator — Monte  Carlo  (HMOMC)  method  to  simulations 
of  underwater  polarized  light  field  under  a  dynamic  ocean  surface.  Field  data  of  wave  slope 
field  and  water  inherent  optical  properties  measured  during  the  RaDyO  Santa  Barbara  Channel 
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(SBC)  experiment  have  been  used  as  input  to  the  HMOMC  calculations.  Simulated  degree  and 
plane  of  polarization  patterns  are  found  to  be  qualitatively  comparable  to  measurements  of 
unpolarized  and  polarized  in-water  light  field  also  made  in  the  SBC  experiment.  Analysis  of 
the  computer  time  for  the  simulation  runs  reveals  that  the  hybrid  model  is  highly  efficient  for 
calculations  of  multiple  frames  of  polarized  light  field  under  a  dynamic  ocean  surface.  A  paper 
on  this  work  has  been  accepted  for  publication  in  Journal  of  Geophysical  Research,  RaDyO 
special  issue  (2011). 

b)  We  have  studied  the  optical  volume  scattering  function  (VSF)  in  a  surf  zone  and  the  inversion 
of  the  VSF  for  derivation  of  sediment  and  bubble  particle  subpopulations.  An  irregular 
hexahedron  model  has  been  used  for  the  calculations  of  VSFs  for  oceanic  particles.  A  paper  on 
this  work  has  been  submitted  to  Journal  of  Geophysical  Research,  RaDyO  special  issue  (2011). 

c)  We  have  investigated  the  effect  of  ice  particle  aggregation  on  light  scattering  by  ice  crystals  in 
regions  of  deep  tropical  convection.  We  developed  an  aggregate  representation  with  10 
ensemble  members,  each  consisting  of  4  to  12  hexagonal  plates.  The  averaged  scattering 
properties  of  aggregate  representations  computed  from  either  the  discrete  dipole  approximation 
or  improved  geometric  optics  method  agree  closely  with  those  from  1000  different  aggregate 
geometries.  This  implies  that  the  aggregation  representation  provides  an  accurate  and  efficient 
way  to  represent  all  aggregates  within  ice  clouds.  A  paper  on  this  work  has  been  published  in 
Applied  Optics,  50,  pp.  1065-1081  (2011). 

d)  We  have  developed  a  new  physical-geometric  optics  hybrid  (PGOH)  method  to  compute  the 
scattering  and  absorption  properties  of  ice  particles.  This  method  includes  analytical 
considerations  of  edge  effects  and  amplitude  and  phase  variations  over  the  geometrical 
wavefront.  Therefore,  it  is  good  for  particles  with  size  parameters  as  small  as  ~20.  A  paper  on 
this  work  has  been  published  in  J.  Quant.  Spectrosc.  Radiat.  Transfer,  1 12,  1492-1508  (2011). 

e)  We  have  compared  three  approximate  methods,  the  surface-integral  method,  the  volume- 
integral  method,  and  the  diffraction  plus  reflection  pattern  from  ray  optics  method,  that 
compute  scattering  of  light  by  dielectric  faceted  particles.  A  paper  on  this  work  has  been 
published  in  J.  Quant.  Spectrosc.  Radiat.  Transfer,  112,  163-173  (2011). 

f)  We  have  made  some  of  the  first  comprehensive  comparisons  between  simulated  and  measured 
underwater  polarized  light  field,  including  both  degree  and  plane  of  polarization,  in  clear 
oceanic  waters.  Satisfactory  consistencies  are  observed  in  the  comparisons.  Based  on  model 
simulation  results,  we  have  also  shown  the  limitations  of  the  in-water  single  scattering 
approximation  in  predicting  the  overall  underwater  patterns  of  degree  and  plane  of  polarization. 
A  paper  on  this  work  has  been  published  in  Applied  Optics,  50,  4873-4893  (2011). 

g)  We  have  made  radiative  transfer  model  simulations  as  a  confirmation  of  the  locations  of 
underwater  neutral  points.  Simulation  results  suggested  that  the  underwater  neutral  points  may 
move  out  of  the  principal  plane  at  longer  wavelength  in  the  visible  when  the  diffuse  skylight  is 
strong  enough.  This  prediction  is  consistent  with  field  measurements.  A  paper  on  this  work 
has  been  published  on  Optics  Express,  19,  5942-5952  (2011). 

h)  We  have  studied  the  effects  of  the  polarization  state  of  an  incident  quasi-monochromatic 
parallel  beam  of  radiation  and  the  orientation  of  a  hexagonal  ice  particle  with  respect  to  the 
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incident  direction  on  the  extinction  process.  It  is  found  that  the  attenuation  of  a  quasi- 
monochromatic  radiation  beam  by  an  ice  cloud  depends  on  the  polarization  state  of  the  beam  if 
the  ice  crystals  are  not  randomly  oriented.  A  paper  on  this  work  has  been  published  in  J.  Quant. 
Spectrosc.  Radiat.  Transfer,  112,  2035-2039  (2011). 

RESULTS 

a)  We  have  used  high  resolution  wave  slope  measurements  and  water  inherent  optical  properties 
(IOPs)  measurements  as  input  to  our  hybrid  radiative  transfer  model,  in  an  effort  to  reconstruct 
the  highly  dynamic  polarized  radiance  field  in  shallow  waters.  Figure  2  shows  snapshots  of 
wave  elevation  fields  reconstructed  from  measured  wave  slope  components  using  a  Fourier 
Transfonn  approach  (wave  slopes  from  Chris  Zappa,  reconstruction  algorithm  from  Howard 
Schultz).  The  wave  slope  measurements  were  made  in  the  Santa  Barbara  Channel  experiment, 
with  a  spatial  resolution  of  roughly  1.3  mm  x  1.7  mm  and  a  frame  rate  of  60  Hz.  However, 
these  high  resolution  measurements  cover  only  a  1  m  x  1  m  surface  patch,  which  is  not 
adequate  for  purposes  of  radiative  transfer  modeling.  We  had  to  extend  this  surface  patch  to  a 
3  m  x  3  m  one  using  periodic  boundary  conditions,  such  that  a  detector  depth  of  1  meter  is 
possible  (i.e.,  the  footprint  of  the  Snell  cone  -  the  region  that  the  sky  radiance  can  reach  - 
viewed  at  that  depth  is  within  the  computational  domain). 


Fig.  2.  Snapshots  of  a  dynamic  wave  elevation  field  reconstructed  from  wave  slope  component 
measurements  in  the  Santa  Barbara  Channel:  (a)  wind  speed  5.1  m/s,  (b)  wind  speed  9.3  m/s.  The 
spatial  resolution  is  roughly  1.3  mm  x  1.7  mm,  and  the  frame  rate  is  60  Hz.  Shown  here  are  1  m  xl 
m  wave  patches.  They  are  extended  to  3  m  x  3  m  and  used  as  input  to  the  hybrid 

radiative  transfer  model. 


Furthermore,  we  have  obtained  relevant  ocean  IOPs  from  RaDyO  collaborators.  Field  data  used  in  the 
model  simulation  include  ocean  particulate  extinction  (cw),  single  scattering  albedo  (  ojp)  and 
backscatter  factor  ( bback )•  The  aerosol  optical  thickness  (AOT)  is  also  used.  A  summary  of  properties 
in  the  simulated  atmosphere-surface-ocean  system  can  be  found  in  Table  1. 
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Table  1.  Summary  of properties  in  the  simulated  ASO  system.  Entries  with  an  asteroid  use  field 
data  from  RaDyO  Santa  Barbara  Channel  experiment.  A  wavelength  of 532  nm  was  assumed. 


Light  source 

Solar  irradiance 

Top-of-atmosphere  solar  radiation  spectrum 

Atmosphere 

Molecular  scattering 
Optical  thickness 

Single  scattering  albedo 
Aerosol  scattering 

Optical  thickness* 

Rayleigh  phase  matrix 

0.112 

1 

Oceanic  aerosols 

0.15 

Air-sea  interface 

Wave  slopes* 

Field  data  under  wind  speed  5.1  m/s  (extended  to  3  m  x  3  m) 

Ocean 

Molecular  scattering 
Extinction  and  albedo 
Hydrosol  scattering* 

Extinction  and  albedo* 
Physical  depth 

Rayleigh  phase  matrix  with  depolarization  ratio  0.047 
cw  =  0.046  m"1,  a>,v  =  0.0457 

Fournier-Forand  phase  function,  where  bhack  is  from  field  data,  combined 
with  reduced  Rayleigh  phase  matrix  with  depolarization  ratio  0.1 
cp  =  0.6522  m'1,  a>p  =  0.9314 

10  m 

Bottom 

Lambertian  surface 

Surface  albedo  A  =  0 

To  resolve  the  fine  angular  distributions  when  real  wave  slopes  are  present,  we  further  improved  the 
hybrid  model  by  treating  the  direct  light  part  (sun  light  getting  into  the  ocean  without  atmospheric 
scattering,  and  then  multiply-scattered  in  the  ocean)  separately.  A  rather  fine  angular  resolution  (1°  x 
2°  in  the  results  shown  below)  can  be  achieved  in  this  part.  The  low-resolution  diffuse  part  (sun  light 
multiple-scattering  in  the  atmosphere  and  enters  the  ocean)  is  then  interpolated  and  added  to  the  direct 
part.  The  total  radiance  field  will  have  the  same  angular  resolution  as  the  direct  part,  which  is 
necessary  to  resolve  the  fine  structures  in  the  angular  distribution  of  dynamic  radiance  fields. 

In  the  simulations,  a  virtual  detector  was  put  at  two  depths,  0.2  meters  and  1  meter.  Figure  3  shows  the 
simulated  angular  distribution  of  downwelling  radiance  field  viewed  at  0.2  meters  (left)  and  1  meter 
(center)  below  the  low-sea  surface  and  at  1  meter  below  the  high-sea  surface  (right).  As  expected,  the 
structures  in  the  angular  distribution  due  to  the  wavy  surface  are  even  finer  when  the  detector  moves 
further  away  from  the  surface  and  when  the  surface  wind  speed  increases.  One  thing  to  notice  here  is 
that  the  radiance  pattern  outside  of  the  Snell’s  cone  shows  almost  no  structure  due  to  extremely  low 
signal. 


Fig.  3.  Simulated  underwater  downwelling  angular  distribution  of  radiance:  Left,  at  0.2  meters 
below  the  low-sea  surface;  center,  at  1  meter  below  the  low-sea  surface;  and  right,  at  1  meter  below 

the  high-sea  surface. 
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Figure  4  shows  the  simulated  degree  of  polarization  (DoP,  left  panel)  and  plane  of  polarization  (PoP, 
right  panel)  viewed  at  0.2  meters  below  the  low-sea  surface.  Unlike  the  radiance  pattern,  these 
patterns  (especially  the  PoP  pattern)  are  highly  dynamic  outside  the  Snell  cone,  since  the  polarization 
states  are  largely  determined  by  multiple  scattering  in  the  ocean.  In  the  DoP  pattern,  wave-like 
structures  are  found  outside  of  the  Snell  cone.  This  implies  that  the  polarization  signal  contains  much 
more  information  about  the  instantaneous  wave  slope  field  than  the  radiance  does. 


1+90 

+45 


0 


I “0.2 
0.0 


Fig.  4.  Simulated  downwelling  angular  distribution  of  degree  of polarization  (left)  and  plane  of 
polarization  (right)  at  0.2  meters  below  the  low-sea  surface. 

In  Fig.  5  we  show  comparisons  between  measured  and  model  simulated  underwater  angular 
distributions  of  radiance,  degree  of  polarization  and  plane  of  polarization.  The  measurements  were 
made  by  RaDyO  researcher  Ken  Voss  in  the  Santa  Barbara  Channel  experiment,  with  a  solar  position 
at  about  55°  zenith,  a  very  clear  sky,  a  surface  wind  speed  of  about  5  m/s,  and  a  detector  depth  of  1 
meter.  In  the  simulation,  we  used  parameters  as  close  as  we  could  get,  and  averaged  the  resultant 
radiance  field  over  every  0.1  seconds  (6  frames)  to  match  up  with  the  integration  time  associated  with 
the  measurements.  Although  it  is  unlikely  that  the  instantaneous  wave  patch  is  the  same  (as  one  can 
tell  from  the  boundary  of  the  Snell  cone),  the  overall  consistencies  in  all  three  patterns  are  surprisingly 
encouraging.  The  model  simulations  failed  to  reproduce  some  of  the  noise  in  the  measured  DoP  and 
PoP  patterns,  which  are  likely  from  both  instrumentation  noise  and  wave  disturbances.  This  could  be 
partly  attributed  to  the  neglect  of  wave  elevations  in  the  model. 
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Fig.  5.  Comparisons  of  DPOL-measured  (upper  panel,  from  Ken  Voss)  and  model  simulated 
(lower  panel)  downwelling  angular  distributions  of  radiance,  degree  of polarization,  and  plane 
of  polarization  at  1  meter.  Sun  is  at  a  zenith  angle  of  about  55°;  the  sky  is  clear;  and  the 

surface  wave  speed  is  about  5  m/s. 


Simulated  underwater  radiances  (unpolarized  light  field)  below  the  high-sea  surface  were  also 
compared  with  its  counterpart  measured  by  RADCAM.  During  this  measurement,  the  solar  zenith 
angle  was  about  45°,  the  sky  was  clear,  and  the  surface  wind  speed  was  about  1 1  m/s.  The  RADCAM 
instrument  operates  at  the  wavelength  of  555  (±20)  mn.  Here  we  are  comparing  against  the  results  at  1 
m.  Figure  6  shows  the  comparison  between  the  measured  (panel  (a))  and  simulated  (panel  (b)) 
underwater  downwelling  radiance  fields  in  units  of  W  m'  mu'  sr'  .  The  measured  radiance  is  slightly 
higher  than  the  simulation  (a  larger  high  radiance  region,  and  brighter  in  the  Snell’s  cone).  Despite 
this  discrepancy,  the  simulation  did  give  dynamic  features  (such  as  distortions  of  the  Snell’s  cone  and 
fine  structures  at  the  edge  of  the  Snell’s  cone)  that  are  qualitatively  similar  to  those  seen  in  the 
measurement.  The  measured  light  field  features  far  more  pronounced  fine  structures  and  a  more 
distorted  Snell’s  cone.  The  sources  of  this  discrepancy  are  two-fold:  1.  the  instantaneous  wave  slope 
field  above  the  measured  light  field  is  different  from  the  wave  slope  measurement  used  in  the  model 
simulation;  2.  the  model  neglected  wave  elevations. 
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Fig.  6.  Comparison  of  downwelling  radiance  field  at  1  m  measured  by  RADCAM  (panel  (a))  and 
simulated  by  the  hybrid  model  (panel  (b))  in  units  of  W  m2  nm1  sr1.  The  solar  zenith  angle  is  45°, 

and  the  surface  wind  speed  is  about  10  m/s. 


b)  In  a  collaboration  with  researchers  under  another  MURI  project,  we  compared  model  simulated 
underwater  light  polarization  with  measurements  made  by  a  multispectral  polarimeter  (CCNY) 
and  by  a  video  polarimeter  (UT).  In  this  study  we  compared  not  only  the  degree  of  polarization 
but  also  the  plane  of  polarization,  both  of  which  play  a  critical  role  for  marine  animals  to  use 
polarization  signals  for  navigation  and  orientation.  Field  measurements  of  water  inherent 
optical  properties,  particulate  scattering  volume  scattering  functions  and  polarized  scattering 
functions  (S12/S1 1  element)  were  used  in  the  simulations. 

Figures  7  and  8  show  the  comparisons  of  the  degree  of  polarization  (DoP)  angle  of  polarization  (AoP) 
at  nine  wavelengths.  The  model  simulations  are  consistent  with  the  measurements  at  shorter 
wavelengths.  In  the  red  the  measurements  are  too  noisy  for  a  good  comparison. 


11 


412  nm 


B1;  SE  =  17.8;  VA  =  0;  z  =  1  m 
440  nm 


(e) 

\  t 

I  J  i\ 

\/+ 

676  nm 


Instrument  Angle  [degree] 


488  nm 


(c) 

555  nm 

(f) 

A 

\  * 

+  jf 

Nt  t 

715  nm 

(i) 

i  i . 

y  l 

I!W|! 

-90  -45  0  45  90  135 


Fig.  7.  Degree  of polarization  (DoP)  patterns  from  CCNY  measurements  and  from  simulations  at 
nine  wavelengths;  the  vertical  dashed  line  indicates  the  angle  that  corresponds  to  a  94°  scattering 
angle;  the  vertical  dotted  line  indicates  the  edge  of  the  Snell’s  window,  with  angles  within  the 
window  to  its  left.  “ SE ”  stands  for  “ Solar  Elevation  ”,  “VA  ”  stands  for  “Viewing  Azimuth  ”,  and  z.  is 

the  detector  depth. 
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Fig.  8.  Same  as  Fig.  7,  but  shows  patterns  of  the  angle  of polarization  (AoP). 
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Figures  9  and  10  shows  ta  comparison  of  model  simulations  against  measurements  from  the  UT  video 
polarimeter  at  the  wavelength  of  5 10  mn.  The  consistency  here  is  not  as  good  since  a  great  deal  of 
estimations  were  made  in  detennining  the  viewing  angles  in  the  measurement.  Nevertheless,  the 
general  agreement  is  encouraging. 


Fig.  9.  Comparison  of  the  degree  of polarization  (DoP)  from  model  simulations  and  from 

UT  video  polarimeter  measurements 
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Fig.  10.  Same  as  Fig.  9,  but  shows  the  comparison  of  the  angle  of polarizatio  (AoP). 
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IMPACT/APPLICATION 


The  HMOMC  code  and  the  fast  irradiance  code  will  become  powerful  tools  to  investigate  the  effects  of 
a  dynamic  wave  profile  on  the  radiance  field,  especially  on  the  polarization  states  and  fast  fluctuations 
of  the  light  field.  The  fast  irradiance  code  will  also  be  helpful  in  understanding  the  fast-varying 
environment  surrounding  ocean  organisms,  which  is  part  of  an  ONR  funded  MURI  project  on 
underwater  camouflage. 

TRANSITIONS 

Due  to  the  efficiency  and  versatility  of  this  new  code,  it  will  be  directly  applicable  to  understanding 
one  of  the  most  formidable  problems  in  global  warming,  i.e.  the  effect  of  broken  clouds  on  the 
reflectivity  of  the  atmosphere.  My  colleague,  Dr.  Ping  Yang  in  the  Department  of  Atmospheric 
Sciences  at  TAMU,  will  use  it  to  interpret  the  measurements  of  the  satellites,  Cloud-Aerosol  Lidar  and 
Infrared  Pathfinder  Satellite  Observation  (CALIPSO),  Moderate  Resolution  Imaging 
Spectroradiometer  (MODIS),  and  Polarization  and  Directionality  of  the  Earth's  Reflectances 
(POLDER).  This  code  can  also  be  used  in  biomedical  studies,  such  as  numerical  simulations  of  the 
light  propagation  in  skin  tissue. 

RELATED  PROJECTS 

We  use  the  results  from  our  other  ONR  Grant  to  use  as  input  to  our  codes  in  our  MURI  study  of 
cephalopods. 
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